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Abstract 

In this work we present a deduction of the Saint-Venant-Exner model through 
an asymptotic analysis of the Navier-Stokes equations. A multi-scale analysis is 
performed in order to take into account that the velocity of the sediment layer is 
smaller than the one of the fluid layer. This leads us to consider a shallow water 
type system for the fluid layer and a lubrication Reynolds equation for the sediment 
one. This deduction provides some improvements with respect to the classical Saint- 
Venant-Exner model: (i) the deduced model has an associated energy. Moreover, 
it allows us to explain why classical models do not have an associated energy and 
how to modify them in order to recover a model with this property, (ii) The model 
incorporates naturally a necessary modification that must be taken into account in 
order to be applied to arbitrarily sloping beds. Furthermore, we show that this 
modification is different of the ones considered classically, and that it coincides with 
a classical one only if the solution has a constant free surface, (iii) The deduced 
solid transport discharge naturally depends on the thickness of the moving sediment 
layer, what allows to ensure sediment mass conservation. Moreover, we include a 
simplified version of the model for the case of quasi-stationary regimes. Some of these 
simplified models correspond to the generalization of classical ones such as Meyer- 
Peter&Miiller and Ashida-Michiue models. Three numerical tests are presented to 
study the evolution of a dune for several definition of the repose angle, to see the 
influence of the proposed definition of the effective shear stress in comparison with 
the classical one, and by comparing with experimental data. 
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1 Introduction 

The Saint-Venant-Exner system (see [17]) is generally used to model the bedload transport 
in rivers, lakes and coastal areas. Sediment transport is usually divided into three types: 
surface creep, saltation and suspension. Surface creep is defined as the type of transport 
where sediment grains roll or slide along the bed. Saltation transport is defined as the 
type of transport where single grains jump over the bed a length proportional to their 
diameter, losing for instants the contact with the soil. Sediment is suspended when the 
flux is intense enough so that the sediment grains reach height over the bed. There is not 
a clear distinction between surface creep and saltation, so that these types of transport are 
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usually called bedload transport. At low Froude numbers, the bedload is the dominating 
transport mechanism which is the regime under study in this paper. 

The Saint-Venant-Exner system (SVE in what follows) is defined in terms of a hy- 
drodynarnical component coupled with a morphodynamical one. The hydrodynamical 
component in most cases is modeled by Saint-Venant system. The equation that describe 
the morphodynamical component is the well known Exner equation, that is a continuity 
equation. Such system can be written under the following form: 

( d t h\ + div x <?i = 0, 


< 


d t qi + div x 


( gi 0 gi \ 

V hi J 



+ gh\\7 x (h -2 + b) + t / p\ — 0, 


(1) 


{ d t h 2 + di v x q b = 0, 

where x — ( Xi,x- 2 ) is the horizontal spacial coordinates, t represents the time variable, 
qi — hi(x,t)ui(x,t) represents the water discharge, h\(x,t) the height of the fluid and 
U\ = (^i[i]A i i[ 2 ]) its horizontal velocity, r is the shear stress at the bottom and pi the 
density of the fluid. The unknown function h 2 — h 2 (x,t) is the thickness of the sediment 
layer (see Figure 1) and qb denotes the solid transport discharge, g is the gravity constant 
and b is the fixed bottom, usually called the bedrock layer. In what follows we shall denote 
g — b + h. 2l being z = p(x, t ) the sediment bed surface. 

To close the system, it is necessary to define the solid transport discharge qb . Several 
formulae for qb can be found in the literature. For example the classical formula proposed 
by Grass [29] assumes that the movement of the sediment begins at the same time as for 
the fluid and both move in the same direction. It is defined by qb — A s |ui| m9_1 Ui where A g 
is a constant which takes into account the grain size and the kinematic viscosity and m g is a 
positive real number, such that 1 < m g < 4. Nevertheless, for practical applications, some 
other type of formulae has been proposed in the literature, for instance, by Meyer-Peter 
& Muller [39], Van Rijn’s [53], Einstein [16], Nielsen [43], Fernandez-Luque & Van Beek 
[6, 20], Ashida & Michiue [1], Engelund & Fredsoe [18], Kalinske [32] or Charru [10]. Such 
formulae are usually presented in nondimensional form and can be written as follows, 

| = sgn(r) (I ^ ) (6 - k 2 9 C )™ 2 (s/0 - h ^9^ , (2) 

where Q represents the characteristic discharge, Q = d s ^/g{l/r — 1 )d s , r is the density 
ratio, r = P 1 /P 2 , being p 2 the density of the sediment particles and d s the mean diameter 
of the sediment particles, ip is the averaged porosity. 

In classical SVE models the sign of qb coincides with the sign of r, the shear stress at 
the bottom. It is usually defined as r = pighiSf, being Sf the friction term. S{ can be 
set by different empirical laws such as the Darcy-Weisbach ( S{ — fui[ui\/8ghi, where / is 
the Darcy-Weisbach coefficient) or Manning formulae (Sf = n 2 ui[ui\/h^ 3 , where n is the 
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Manning coefficient). In general we can write r = p\g((hi)u\\ui\, being ((hi) a function 
depending on the considered friction law. 

The Shields stress, 9, represents the ratio between the agitating and the stabilizing 
forces on a sediment grain in the bed, 


e 


V\d 2 s 

g(p 2 - pi)dl' 


(3) 


and 9 C is the critical Shields stress for incipient motion. 

The constants ki, mi, l — 1,2,3 are positive real numbers depending on the model. 
Being usually at least one of the parameters mi, m 2 or m 3 equals to zero. For example, 
Meyer-Peter & Muller’s model is defined by 


and Ashida & Michiue’s model is defined by 

^ = sgn(r) (9 - 9 C )+(V9 - sj¥ c ). (5) 

Finally, by (•) + we denote the positive part. In equation (2) the positive part implies 
that the sediment moves when the modulus of the shear stress is bigger than a given critical 
value. 

Although the classical SVE model is largely used, it presents several disadvantages: 

(i) The SVE model has not a dissipative energy equation associated to the system. 

(ii) Solid transport discharge formulae are derived by using the hypothesis of nearly 
horizontal sediment beds, that is, V x r] 0. 

(iii) Solid transport flux is independent of the thickness of the sediment layer. Thus, the 
mass conservation property for the sediment given by third equation in ( 1 ) may fail 
(see [40]). 

Concerning the first item, we should remark that there exist in the literature some 
simpler solid transport formulae for which the corresponding SVE model has an associated 
dissipative energy equation (see for example [36]). Nevertheless, up to our knowledge, no 
general result exists in the bibliography in this sense. 

The second item implies that classical formulae cannot be used in several problems of 
interest (see [33]) because they are derived by using the hypothesis of nearly horizontal 
sediment beds. 

As mentioned before, the Shields parameter is the coefficient between agitating forces 
and the stabilizing forces. Classical formulae consider that the only agitating force is the 
bottom shear stress, concretely |t| d 2 s . Nevertheless, in the experiments presented by Lysne 
in [37] it can be seen that gravity is another contributing factor as an agitating force (see 
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also [27]) for sloped sediment beds. Then, it is necessary to take into account gravitational 
forces in order to obtain a solid transport discharge that can be applied in arbitrarily 
sloping beds. 

This has been done in the literature in several ways. For instance, the simplest way to 
take into account the sediment bed slope in the definition of the solid transport discharge 
is to include a diffusion term. Engelund and Hansen proposed in [19] a formula that can 
be written under the following form, 


q b = fc|ui| m - cV x r^ , (6) 

k, m, and c being constant parameters of the model (see also [58] and [59]). Equation (6) 
can be seen as a modification of the Grass model. An adaptation of this formula for curved 
channels was proposed by Struiksma et al. in [57]: 


q b = k\ui\ m (1 


cd s (b + h 2 )) sgn 




s being the streamwise coordinate and f s the shape factor of the grains (see [57] for more 
details). Note that this definition implies that the direction of the sediment transport does 
not coincides with the direction of the velocity of the fluid. The direction is determinated 
by the sign of the vector (u\ — j-gS7 x ri). Note that 


sgn 



Ul ~ -jrgVxV 
ui - jjg'Vx'n 


cos a 
sin a 


a being the angle of the transport direction, where 


u i [2] r e d X2 ri 

tana =-G—-—. 

U 1 [ 1 ] — f s g Ox\ V 


(7) 


Let us remark that in [57] the direction of the sediment transport is defined in terms of 
its angle, given by (7), rather than defining it in terms of the sign vector. Authors also 
include a correction of the angle due to the transversal velocity in curves. For the sake of 
brevity we do not include a discussion on the modification of the transport angle in curved 
channels, since it is not the aim of this paper. 

A more used extension of classical formulae to arbitrarily sloping sediment beds is to 
consider a modification of the critical Shields parameter, replacing 6 C by 6 C (see [27] and 
[20]). For the scalar case, ID flows, this modification can be written as follows: 


9 


C 


0c 



sgn(r) 
tan 5 



9 c + d sgn(r) d x ri, where •& = — 

tan o 


( 8 ) 


In the works of Kovacs and Parker [33], Seminara et al. [56] and Parker et al. [47] several 
extensions for the vectorial case are presented, where the computation of 9 C takes into 
account lateral slopes. 


5 



Note that usually in the definition of the modified critical Shields parameter we do not 
find the sign of r. This is due to the fact that formulae are usually presented for the case 
of positive velocities only. 

In previous definition tan S is the friction coefficient corresponding to the internal fric¬ 
tion angle of the material. Typical values of i) and 9 C are, $ — 0.1 and 9 C — 0.047, what 
implies that 9 c /i) = tan 6, with 8 25°, as proposed by Fredsoe in [27]. The angle 8 = 25° 

is lower than the repose angle close to 32°, although lower values have also been suggested 
(see [10] and references therein). 

Let us remark that the definition of 9 C in (8) is based in two arguments: first, grav¬ 
itational forces due to the sediment bed slope are incorporated as agitating forces in the 
definition of the effective Shields parameter, 


9 P ff = 


\r e s\d 2 s 


g(p2 ~ Pi)d 3 s 


with T eff = r - $(p 2 - pi)gd s d x (b + h 2 ), 


(9) 


where r e g is the effective shear stress. Second, rather than replacing 9 by 9 e g, it is usually 
assumed that this is equivalent to replacing 9 C by 9 C , defined by (8) (see [27]). More 
explicitly, it is usually assumed that (9 e g — 9 C ) — (9 — 9 C ). 

Nevertheless, this second assumption is not true in general, 


9 e & ~ 9 C ^ 


\r\d; - #sgn(r)g(p 2 - pi)d 3 d x {b + h 2 ) 
g{p2 - Pi)d 3 s 


— 9 C — 9 — 9 r 


The problem arises from the fact that the absolute value (or the norm in the vectorial case) 
is neglected in the definition of 9 e g (see for example [10] and [27]), which should be taken 
into account, depending on the sign of r and r e g. 

In fact, Fowler et al. proposed in [26] a modification of the Meyer-Peter & Muller 
formula that consist in replacing 9 by 9 e g, instead of replacing 9 C by 9 C . The model proposed 
by Fowler et al. can be written as follows: 


Ob 

Q 


8 sgn(r e ff) 


ha 

ha 




3/2 

+ 


where ho is an averaged value of the thickness of the sediment layer. The fact of introducing 
an explicit dependence on h 2 in the formula is also interesting. Indeed, this allows to ensure 
the mass conservation property for the sediment layer (see also [40]), which was the second 
issue noted previously as a disadvantage for classical SVE models. 

Related to this problem, in [23], Fernandez-Nieto et al. introduce a modified general 
definition of the solid transport discharge for SVE models that takes into account the 
thickness of the sediment layer. Then mass conservation is ensured. Moreover, the proposed 
formula has the advantage that it reduces to a classical solid transport discharge formula 
in the case of quasi-uniform regimes. This is in fact the regime where usually classical 
formulae are derived. 
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The main objective of this paper is to show that the SVE model can be deduced through 
an asymptotic analysis of the Navier-Stokes equations. Moreover, we obtain that following 
this process some improvements on the classical SVE model are reached: (i) We obtain 
that the deduced model verifies exactly a dissipative energy equation. This deduction 
and the proof of energy allows us to understand why classical models do not have an 
associated energy. Moreover, it allows to prove that the classical SVE model can have an 
associated energy by introducing a simple modification, (ii) In the deduction process of 
the model the pressure terms introduce a modification that must be taken into account 
if we consider applications where the sediment bed is not nearly horizontal. We also see 
that this modification coincides with some of the alternatives proposed in the literature for 
some special cases, (iii) The solid transport flux depends on the thickness of the moving 
sediment layer, then mass conservation is ensured. 

The paper is organized as follows: In Section 2 we present the 3D system of equations 
considered as starting point from which the models are deduced, the proposed models and 
the results on the associated energy. The formal deduction of the models by an asymptotic 
analysis from the Navier-Stokes equations is detailed in Appendix A. This is done by 
developing a multi-scale analysis in space and time. Section 3 is devoted to numerical 
tests. In the first test we study the evolution of a dune for different values of the the 
repose angle. The purpose of the second test is to compare the influence of the deduced 
modification in the definition of the effective shear stress respect to classical ones. A 
comparison with experimental data is presented in the third test. Finally, the conclusions 
are presented in Section 4. 


2 The asymptotic Saint-Venant-Exner models 

Following an asymptotic analysis, we derive two main models which are presented here (see 
Appendix A for details). The difference between these models is the friction law used at 
the fluid and sediment interaction level. In Subsection 2.1 we summarize the starting 3D 
system of equations and the hypothesis considered for the derivation of the models. The 
models deduced in this work are shown in Subsection 2.2. First, they are presented with 
the same notation as they are deduced (see Appendix A). Secondly, they are rewritten in 
terms of the Shields parameter or the effective Shields parameter. The associated energy 
for these models is presented in Subsection 2.3. Moreover, we include in Subsection 2.4 
a result that justifies that classical SVE model may have an associated dissipative energy 
provided a second order correction is made in the friction term appearing in the momentum 
equation. 

2.1 The 3D initial system 

We consider two immiscible layers of different materials with different physical properties: 
velocity, pressure, density and viscosity. The two layers are related through the interaction 
terms at the internal interface levels. In the following subsections the starting systems 
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of equations, the definition of the physical domain and the considered hypothesis for the 
deduction of the models proposed in this paper are specified. 

2.1.1 Physical domain and governing equations 

We consider a cartesian coordinate system where x represents the horizontal 2D direction 
and z the vertical one. Let us define the physical domain for the fluid and sediment layers 
by Of (t) and D 2 (t) respectively; t being the time variable. Usually in the context of bedload 
transport it is assumed that the sediment domain can be decomposed into two layers: one 
that moves due to the action of the convection of the upper fluid, the mobile sediment 
layer with thickness h rn , and a second one composed of sediment that is not moving but is 
susceptible to come into motion, with thickness hf. This leads us to define four boundaries 
in the domain (see Figure 1): 

r s = {(x,z) G M 3 /x Euj, z = rj s (x, t)}( 

Tig = {(x,z) G M 3 /x G ui, z — rj(x,t)}-, 

Tf — {(x, z) G M 3 /x G uj, z — r]f(x , £)}; 
r& = {(x, z) G M 3 /x G uj, z = b(x)}. 

where w is a domain in M 2 . The water free surface is defined by z — rj s (x,t), where 
r/ s (x,t ) = hi(x,t) + h 2 (x,t) + b[x) and the fluid/sediment interface by z — r/(x,t), where 
rj(x, t ) = b(x)+h 2 (x, t). hi(x, t ) denotes the height of the water column. The sediment layer 
is decomposed as h 2 {x,t) — h m (x,t) + hf(x,t ), and then the internal sediment interface is 
z — r/f(x,t), where r)f(x,t ) = b(x) + hf(x,t ) (see Figure 1). 

Thus, we consider a time-dependant domain Q(f) = Di(t) U D 2 (t) U T^ U ri j2 (t) U T s (t), 
being: 

Di (t) = {(XfZ) G M 3 /x G uj, r](x,t ) < z < r] s (x,t)}; 

n 2 (t) = n 2j /(t) u n 2)m (t); 

where 

fl 2 j(t) — {(x,z) G M 3 /x G uj, b(x) < z < r]f(x,t)}-, 

D 2 , m (t) = {{x,z) G M 3 /x G uj, rjf(x,t ) < z < rj(x,t)}; 

Moreover, let us denote by T m the mass transference between the static and the mobile 
sediment domains, floj and D 2)TO . T m is defined as the difference between the erosion rate 
(i e ) and the deposition rate (id), T m = z e — id (see [23] and Section 2.2.3). 

As a general rule in notation, we will use the subscript 1 to denote the upper layer 
(fluid) and the subscript 2 for the lower layer (sediment). We denote by 

Vi = (u h Wi) 

the velocity field for each layer with Ui — ( Ui p], u % [ 2 ]) the 2D horizontal velocity. We denote 
by pi the density and p t the pressure. Moreover, //* and = fa/pi, denote the dynamic 
and kinematic viscosity coefficients respectively, for * = 1,2. 



2 = V. 




z = q 


\ z = Vf 


Figure 1: Sketch of the domain for the fluid-sediment problem 


For each layer (* = 1, 2), we start from the 3D Navier-Stokes equations for incompressible 
fluid and sediment components: 



( 10 ) 


In this system g represents the gravitational vector and a l the stress tensor associated to 
each layer. 

To complete this system, we must give the stress tensor expressions, the interactions at 
the internal interface levels (Ti i2 and Tf), as soon as boundary and kinematic conditions. 
They are specified in the following subsection. 

2.1.2 Closures 
Stress tensors 

We shall define the stress tensors as follows: 


Oi = a\ — pjId, for i — 1,2; 


where a[ is the deviatoric part of the stress tensor. 

We assume that the stress tensor for the fluid layer follows a Newtonian form with constant 
dynamic viscosity and then a\ is given by: 


cr) = 2 j U 1 L>('Ui), 


D(v) being the rate of deformation tensor, D(v) — |(Vu + V t v). 

For the sediment layer we consider a non-Newtonian rheology. Recent works have been 


devoted to demonstrate through experimental results the resemblance between the bed 
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load transport phenomena and the granular flows behavior [45, 46, 2, 52], Ouriemi et 
al. proposed in [45] a two-phase model for bed load transport in laminar flows. Where 
a Newtonian law for the fluid phase and a frictional rheology for the particulate phase, 
namely a Coulomb type friction is considered. Chauchat et al. (see [7], [ 8 ]) study a three- 
dimensional two-phase model where a Drucker-Prager rheology describe the stress tensor 
associated to the particulate phase. This definition is consistent with the Coulomb friction 
model proposed in [45]. A comparison with experimental data for bed load transport can 
be seen in [ 2 ]. 

Following these references we define 

<y'i = 2/x 2 D(v 2 ), 

with a non-constant viscosity /j 2 , given by a Drucker-Prager model. By using a regular¬ 
ization of the Drucker-Prager model (see [35]), we can assume that /j 2 is expressed as a 
function of D(v 2 ) and p 2 . 

Nevertheless, let us remark that in order to deduce a first order Reynolds-type model 
for the sediment layer, we do not need further specifications about this rheology. In fact, 
the depth-averaged model will be written in terms of the boundary conditions, that is, the 
Coulomb friction law, which is compatible with the Drucker-Prager rheology. 

Friction laws 

We must define the friction laws at the mobile-static sediment interface and at the fluid- 
sediment interface. 

o Friction law at the mobile-static sediment interface 

As explained before, a Coulomb friction condition at the interface between the static 
and the moving sediment particles, z — r)f(x, t), is considered. Denoting by Nf the unitary 
normal vector to the interface T/ and by 6 the repose angle, we write the Coulomb friction 
law as 

(a 2 N f ) T = -(sgn(u 2 )tanJ ((cti - a 2 )N f ) • A r /j _ , (11) 

where Nf — (—V x r/f, 1 f + |V.r? 7/| 2 and the subscript (• )p denotes the tangent com¬ 
ponent of a vector. 

Let us remark that the use of a Coulomb friction can also be interpreted as a mechanism 
to approximate the collision effects of saltating grains in the computation of the bedload 
transport formula (see [33]). 

o Friction law at the fluid-sediment interface 

At the fluid-sediment interface the interaction between fluid and sediment is defined 
through a friction force: 

(aiNflr = (<ToNflr = (fric, O) 4 , ( 12 ) 
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where N v = (-V x rj, 1 Y/y/l + |V X ??| 2 is the normal vector at the interface pointing from 
layer 2 to layer 1. 

In particular, we consider two classical friction laws for which we obtain two different 
models that can be written under the same structure. 

• Linear friction law: 

fric = C («i — u 2 )p=7?) (13) 

where the coefficient C has velocity dimension. Moreover, taking into account the 
results presented in [42] (see also Remark A.2), C can be assumed proportional to 
h m . Following the analysis of Seminara et al. [56], the drag coefficient is proportional 
to tai\{5)/9 c . That is, inversely proportional to d, defined in (8). Taking into account 
these remarks we may define 

c = (1/r - 1 )gh m 

?V(1 /r - 1 )gds 



• Quadratic friction law: 


fric = Ci |(«i - u 2 )\ z=v \{u 1 - u 2 )\ z=v . (15) 

In this case C\ must be adimensional, so taking into account previous arguments on 
the definition of the drag coefficient, we can define 


Ci 


d'm 

tfd s ' 


(16) 


Boundary and kinematic conditions 

The following boundary and kinematic conditions are imposed at each interface to complete 
the system: 

• At the free surface, z = ri s (t, x ) = b(x) + h 2 (x, t ) + h\(x, t ): 

— The surface tension condition: (cri • N s ) n — 0 where N s — , 1 =(— V x r] s , 1)* 

^ 1+\V x ris\ 2 

is the unitary outward normal vector to the free surface and the subscript n 
denotes the normal component. 

— The kinematic condition: d t j] s — V\ ■ N s . 

• At the fluid/sediment interface, z — r](t,x) — b(x) + h 2 (x,t): 

— The kinematic conditions corresponding to both velocities: 

d t r) = Vi ■ N v = v 2 ■ N v . 

- The continuity of the normal component of the tensors: (<7i • N v ) n = (a 2 ■ N v ) n . 
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- The friction law defined by (12)-(16). 

• At the internal sediment interface, z = r]f(t,x) — b(x) + hf(x,t ): 
— The conservation of the sediment mass, 


9tVf =v 2 - N f - T m ; 


T m being the mass transference term (see Section 2.2.3). 

— The Coulomb friction law defined by (11). 

• At the bottom, z — b(x): 

— The no penetration condition: v 2 ■ N b = 0, where the unitary normal vector to 
the bottom is N b — (— V x b, 1)*/-y/l + \V x b\ 2 . 

2.2 Proposed models 

In this section we present the final models obtained through an asymptotic analysis of 
the 3D system (10). Following the work performed in [24], we derive a mathematical two 
dimensional SVE type model for bed load transport. Thus, the models are composed by 
three equations: the first two equations represent the Saint-Venant system used for mod¬ 
eling the upper fluid layer and the third one describes the evolution of the moving bed 
through a lubrication Reynolds equation. 

As it is classically considered in such kind of models, we take into account two different 
time scales for the hydrodynamics and the sediment evolution. Therefore, the ratio be¬ 
tween the hydrodynamic and morphodynamic time scales, defined by e 2 , being small. We 
obtain first order models, so we neglect the terms of order e 2 as commonly done to obtain 
SVE type models, (see [28] and references therein). 

The two models obtained depends on the considered friction law at the fluid-sediment in¬ 
terface, given by (13) and (15) respectively. The complete derivation of these models is 
detailed in Appendix A. 

First, we present the models deduced under a unified formulation. Then, we analyze 
each of them in order to relate the sediment velocity with the classical motion threshold in 
terms of the Shields parameter and the effective Shields parameter. Finally, the closure of 
the systems is discussed and simplified models corresponding to a quasi stationary regime 
will be presented. From now on, we use the superscript (LF) to denote the properties of 
the model deduced with a linear friction law, and ( QF ) for the model using a quadratic 
friction law. 

The model deduced in this paper can be written under the form of a SVE system form 
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as follows: 


d t h] + div^gi = 0, 

d t qi + div x (/ii(ui <g> iq)) + ^.gV x /? 2 + gh x W x {b + h 2 ) + = 0, 

d t h 2 + div x (/i m -\/(1/r - l)c/d s ) = 0, 

„ d t hf = -T m . 


with 

P = V x (rhi + h 2 + 6) + (1 - r)sgn(n 2 ) tan 6. 
and the velocity Vb defined for each friction law as 

1 i? 


(LF) 

Vl = 


Vi 1 / 7 ' - l )ad : 


XQF) 


Vi 1 / 7 ' - 1 )gd, 


■■U 1 


= «1 - 1 

1 — r 

|P| 1/2 sgn(P). 


t? y / 2 

1 - r) 


(17) 


(18) 


(19) 

( 20 ) 


We can write these velocities in terms of the Shields parameter. For this aim, we use now 
the expression of the shear stress and the definitions of the friction coefficients C and C\ 
given by (14) and (16). We obtain, 

4 LF) = S go(xr>)(«S F| - <y + , 

where r e g and £/)g F ^ are defined by (24) and (25), respectively, for the case of a linear 
friction law. For the case of a quadratic friction law v^ 11 is defined by (35). If we consider 
a linearization of the quadratic friction law (see Subsection 2.2.2), we obtain 

-r= ^rwrp 2 - <>y\, 

where r e ( g F) and 6^ F) are defined by (32)-(33). 

In the following subsections we describe the deduction of previous definitions of Vb for 
each friction law. 


2.2.1 Model deduced with the linear friction law (13)-(14) 

From the asymptotic approximation at first order, we introduce the shear stress and the 
Shields parameter as 


t (lf) 

Pi 


Cu\ and 


(1/r - 1 )gd s 


C\ Ul | 

(1/r - 1 )gd a ‘ 
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Using the definition of the friction coefficient C in equation (14), we have 


u l 


._= = -dsgn(ni)-^ 0^ LF \ 

\/(l l r - 1 )gd s h m 

Thus, taking into account (18) and the definition of i) in (8), the velocity of the sediment 
layer reads in this case 


v b LF) = ^sgn(ui)^ 1 9 (LF) - —V x (r/ii + rj) - sgn [u 2 )9 c . 
hm 1 — r 


( 21 ) 


o Influence of the Coulomb friction law. 

Note that the sign of the velocity of the sediment layer, sgn(u 2 ), has still to be defined. 
Observe that this coefficient comes from the contribution of the Coulomb friction law at 
the interface between moving and static sediment particles (see (18)). In order to specify 
the sign of U 2 , we remark that Coulomb friction force has the same sign of the net force 
acting on the sediment. That is, 

sgn i _ r )sgn(u 2 )tan<^) = sgn (^==== “ 7 + h 2 + b )) • 

\ i_r / K v0-/ r - l )gd s l ~ r ' 


Then, using that I) — 


sgn(w 2 ) = sgn ( Ul ==== = - —*—'V x (r/ii + h 2 + 6) ] • (22) 

Vv / ( 1 A- 1 )5^ i_r / 

Furthermore, Coulomb friction force implies also a threshold in order to stop the motion 
of the sediment layer when the net forces acting on the sediment is not large enough to 
compensate the friction force. In this case, we obtain that the velocity of the sediment 
layer Vb must be zero under the following condition: 


h m U\ 


\J (l/ r - 1 )gd s 


D h m 
1 — r 


V x (rhi + h 2 + b) 


< 


9 h m 
1 — r 


(1 — r) tan <5 = h m 9 c . (23) 


o Link with the effective and the critical Shields parameter. 

The stop criteria (23) gives us the way to find the relationship with the classical motion 
threshold. Thus, we introduce the modified shear stress including the gravitational forces, 
called the “effective shear stress”, defined in our case as: 


fl LF ) 

T eS 

Pi 


Kd s fl LF ^ 

h rn Pi 


—+ h 2 + b). 


(24) 
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The corresponding definition of the effective Shields parameter is 

« F| - tSt - wSk, {*«**£*” ~ +*»+•>) - 

(25) 

In particular, with these definitions we obtain that sgn^jj ) = sgn^) thanks to (22) 
and the definition of fi LF \ As a consequence, from (21), the velocity of the sediment layer 
can be written as 

v ( b LF) = sgn (t^ f) )( 9^ f) -6 c )+. (26) 

Observe that condition (23) coming from the Coulomb friction law is equivalent to the 
classical one, > 9 C . 

Remark 2.1 (Comparison with the classical effective shear stress) The definition 
of the effective shear stress ( 24 ) can be seen as a generalization of the classical one, given 
in equation (9) that can be written in vectorial form as follows (see [26]), 


— = — - d(l /r - 1 )gd s V x (h 2 + b ). 
Pi Pi 


(27) 


There are two main differences: the first one is the factor (. nd s /h m ) multiplying the shear 
stress, which comes from the fact that the model is expressed in terms of the thickness of 
the moving sediment layer. The second main difference appears in the second term, that is, 
we consider (gd s dV x (rh\ + hi + b)/r) instead of (gd s d(l/r — l)V x (/i 2 + &))• Nevertheless, 
both terms are related: 


gd s d(l/r - l)V x (h 2 + b) = 


gd s d 

r 

gdfid 

r 


V x (rhi + h 2 + b) - rV x (h\ + h 2 + b) 


V (rh\ + h 2 + b) - gd s KV x (hi + h 2 + b). 


(28) 


Then, both definitions coincide in the case of constant free surface, V x (h\ + h 2 + b) =0, 
but could be relevant otherwise. 


2.2.2 Model deduced with the quadratic friction law (15)-(16) 

For the quadratic friction law we have the following definition of the fluid shear stress for 
the first order approximation 


r(QF) = Cfiufux, and then e {QF) = 


(29) 


Pi (1/r — l)gd s ’ 

Then, using the definition of C\ in (16) and again (18) and (8), the velocity of the sediment 
is written as 

i / 2 


vl QF ' > — sgn(ui)fi?^ 


d 


1 — r 


Vxfrhi +h- 2 + b) + sgn (u 2 )9 c 


1/2 


sgn 


d 


1—7’ 


V x (r/ii + h 2 + b) + sgn(« 2 )^< 


(30) 
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Following analogous arguments as in the case of linear friction law, we can define 


sgn(u 2 ) = sgn(<f>) 


( 31 ) 


where 


$ = sgn(ui) 


d 


1 — r 


V x {rhi + h 2 + b) 


1/2 


sgn 




1 — r 


and similarly to (23), the stop criteria is 

\h m $\ < h m 6 c . 


V x {rhi + h- 2 + b)j , 
(32) 


For this model we introduce the effective shear stress as 

rg F) = giX/r - l)d a |$ 1$ and 9«P = |<F| 2 , (33) 

where <t> is dehned by (32). The former stop criteria, derived directly from the model with 
the Coulomb friction law, becomes 

0S F) > 0c (34) 

Observe that in this case we are not able to obtain an explicit expression of the velocity 
v b QF ^ in terms of 0^ F \ due to the presence of the modulus and the sign vectors. Never¬ 
theless, the definition of must take into account previous criteria. Thus, thanks to 
(30) and (34), we give 

(qf) _ f v ( b QF) defined by (30), if 6^ F) > 9 C , 

Vb ~ \ 0 , otherwise. W 

We may also consider the following approximation, which can be seen as a linearization, 
of the quadratic friction law (15): 

fric C\ |(^l ^ 2 )\z=r/\{,^l ^ 2)1 z=7] T f-U |(^ 2 ^ 2 )p =^|(^2 ^ 2)1 z—rjf 

being u* 2 the velocity of the system without considering the Coulomb friction law at the 
internal interface of the sediment layer. In this case, we obtain the following bedload 
velocity: 

c f> =- «; /2 ) + - oo 

This corresponds to defining a solid transport discharge with a structure similar to that of 
classical formulae (see for example the Ashida & Michiue’s model (5)). 

Remark 2.2 Definition (33) of the effective shear stress is different of the classical ones. 
Nevertheless, as in previous case (see equation (24)), the main difference is to consider the 
term V x (r/q + h 2 + £>)/( 1 — r) instead ofV x (h 2 + b). Note also that a linearization of (33) 
leads to previous definition of the effective shear stress ( 24 ). 
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2.2.3 Sediment mass transference and simplified models 

To close the models, it is necessary to define T m , the mass transference between the moving 
and the static sediment layers. It is defined in terms of the difference between the erosion 
rate, z e , and the deposition rate z d , 


z e Zd- 


Following the definitions used in [23], 

i e = K e Q (fl-fl c ) + , z d = K d Q h %, (37) 

d s { 1 - <£>) dj 

I\ e and K d being the erosion and deposition constants respectively. 

Nevertheless, it is necessary to remark that in such definition, the erosion rate does 
not take into account gravitational effects appearing in arbitrarily sloping beds. Thus, for 
the case of sediment beds which are not nearly horizontal, we propose to replace 6 by 9 e g. 
Thus, we set 

Ze — dC e 

with (9 e ff defined by (25) or (33). 

A simplification of the proposed model is to consider a quasi-uniform regime, where 
the deposition rate equals the erosion rate, that is T m — z e — z d = 0. In this case, from the 
definition of z d and i e , we have 

hm= A e eS -9 c ) + - (38) 

K-d\ 1 — V) 

If we introduce this definition of h rn in (17), we obtain the following system 
( d t hi + div^i = 0 , 


(l-^° eS 9c)+ 


I d t q! + div x (hi(ui <g)Ui)) + ^gV x h\ + gh 1 \7 x (b + h 2 ) + = 0, ^ 

d t h 2 + div x ^ (^eff - 0c)+ v b \J (1/r - 1 )gd^j = 0, 

with v b and 9 e g defined by (24)-(26) for (LF) and by (30)-(33) for ( QF ). 

Let us focus on the case of a linearization of the quadratic friction law, for which v b is 
defined by (36). Then, we obtain a SVE model with the following definition of the solid 
transport discharge, 

^ = Sgn(Teff) ^ ^ (fleff - 9 C )+ (V^S ~ V^) • ( 40 ) 
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This definition can be seen as a generalization of the Ashida-Michiue model for arbitrarily 
sediment sloping beds, see equation (5). 

Other possibilities for the definition of h rn can be found in the bibliography. For in¬ 
stance, it is possible to define (see [56]) 

{9 eS - 0 C )T~ 

m (i-vOVS^-WE’ ( j 

k\ and k 2 being two parameters of the model. Note that (41) is an approximation of (38), 
which is valid except near the threshold of motion. 

Let us also remark that for the case k- 2 — 1 and Vb defined by (36) we obtain a SVE 
model with the following definition of the solid transport discharge, 

§=sgn(TjT)^(0..-fyr. (42) 

This definition can be seen as a generalization of the Meyer-Peter&Miillcr model for arbi¬ 
trarily sediment sloping beds, see equation (4). 

Remark 2.3 The simplified model defined by (42) coincides with the classical Meyer- 

PeterSzMiiller model for k\ — 8 and by neglecting all gravitational terms in the definition 

/ \ !/2 

of T e g. That, is, for T e g defined by (33) with $ = sgn(ui)rd^0^ QF " ) J , where is 

defined by (29). 

We obtain also that with this definition of T e g the simplified model (40) coincides with 
the classical Ashida-Michiue model when K e /Kd — 17. 

By analogous arguments we can deduce a large range of classical models. What allows 
us to see these classical models as simplified models which can been deduced from an 
asymptotic expansion of the Navier-Stokes equations, by considering the coupling between 
a shallow water layer and a Reynolds one. Then, for example - see previous remark-, 
we can deduce that Meyer-Peter&Muller model corresponds to define h m by (41) with 
k\ — 8 and k 2 — 1. By while the Ashida-Michiue model corresponds to define h rn by (38) 
with K P JK r ] = 17. In Figure 2 we can see the comparison of both definitions of h rn for 
Q c — 0.047. Where we can observe that both definitions of h m are close. 

2.3 Energy balance of the models 

In this subsection we present the main result regarding the energy balance associated to 
the models ( LF ) and ( QF) presented before. In particular we prove that they admit an 
exactly dissipation energy. The detailed proof of the theorem is given in the Appendix B.l. 

Theorem 2.1 The model (17) has a dissipative energy balance. More explicitly: 
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Figure 2: Definitions of h m corresponding to the Meyer-Peter&Miiller and the Ashida- 
Michiue models. 


for the ( LF ) model: 


dt (\ r 9(hi + h 2 + b) 2 + ^rhi\u\\ 2 + gh 2 x(l - r)sgn(u 2 ) tan <5 
+div x ^'rhiui (yy~ + 9(hi + h 2 + b)j ^ 

+div x {h m u 1 V - hmyyPV \J(l/r - 1 )gd^j 

< -ghmzr^-z \J (i/ r - l )g d s\v\ 2 , 


1 — r 


for the ( QF ) model: 


d t (yrg(hi + h 2 + b) 2 + yh x \ui\ 2 + gh 2 x(l - r)sgn(u 2 ) tan <5 
+div x ^rhiutlyy + g(hi + h 2 + 

+div x (h m uiV - h m \JYyfPl'Pl 11 ' 2 sgn(V) \J (1/r - 1 )gd t 
< Vi 1 / 1 "- l )9ds\V \ 3/2 , 


(43) 


(44) 


where V = g((rh\ + h 2 + b) + x (1 — r)sgn(u 2 ) tan5) and V is defined in (18). 

2.4 Energy balance of the classical Saint-Venant-Exner model 

As mentioned in the introduction there exist in the literature some simple solid transport 
formulae for which the corresponding SVE model has an associated dissipative energy 
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equation. Nevertheless, up to our knowledge, no general result exists in the bibliography 
in this sense. 

In this subsection we present a result that shows that by including a second order 
correction in the definition of the friction law in the momentum equation of the Saint- 
Venant system, we obtain a dissipative energy balance for any classical SVE model whose 
solid transport discharge can be written under the general form ( 2 ). 


Theorem 2.2 Let us consider the general SVE system (l)-(3) and the following definition 
of the friction term in the momentum conservation equation of the Saint-Venant model, 

r/pi = KOiWM + ^€m\/gC(hi)'R~ 

We have that if _ __ 

Vd - \[d c + sgn{r)gV7 X {rh x + h 2 + b) > 0, 

then the SVE model satisfies the dissipative energy balance: 

d t Q r g{h t + h 2 + b ) 2 + ^rhi\ui \ 2 - h 2 sgn(r)\/W c 

+ div x + g(h i + h 2 + b))^j 

+ div x (ZmU!\jgC,(hi)K - f m TZlZ \/(l/r - ljgd. 

< -5C(^l)l«l| 3 - tmVi. 1 / 1 " - l )gds W 2 


where 


£>m 


h e m1 (6 - Me):T 2 (v^ - hVo c y;*- 


(45) 


1^ -v (y/e-y/Bc + s^gVxlrhi + hv + b))- 

7 Z — sgn{T)\fff c — gW x {rh\ + h 2 + b) and 7 Z = sgn{r)x\fT c — g(rhi + h 2 + b). 


The proof of this theorem is similar to the one of Theorem 2.1 and can be found in the 
Appendix B.2. 


Remark 2.4 Note that the considered modification of the classical SVE model introduces 
a second order term in the system. That implies that the classical SVE model verifies a 
dissipative energy equation up to second order terms. 


3 Numerical tests 

In this section we present three numerical tests, for the simplified model obtained from 
(17) with (41), that is, the solid transport equation reduces to (42). This coincides with a 
classical Meyer-Peter&Mullcr model with a modified shear stress given by r e s The purpose 
of these tests is to study the influence of the effective shear stress. The purpose is to study 
what are the effects of the different terms that appear in this modified shear stress. 

The numerical results follow from a combination of the scheme described in [ 6 ] with a 
discrete approximation of bottom and surface derivatives. The numerical simulations are 
done with a CFL number equal to 0.5. 
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3.1 Test 1 


In this first test we propose to study the influence of tand = 0 c /i? in (24). To do so, let us 
consider the following initial condition 

M0,x) = | ^henriae}’ h ^ x ) + MM) = MM) = 1.5. (46) 

The initial condition is shown in Figure 3(a). We set the boundary condition q(t, 0) = 1.5. 
The parameters for the model has been set as follows 

r = 0.34, d s = 10“ 3 , 

<p = 0, K e = 10, K d = 1, e c = 0.047 

with the Manning friction law and as Manning coefficient n — 0.01. The computational 
domain used is [0,10] with 800 points. 

Remark that this is a rather severe test: the bottom is discontinuous and thus we have 
d x {h 2 + b) — ±oo initially. We run the test for different values of 8 ranging from 25° to 89°. 
In Figure 3(b) we show the surface and bottom at time t — 2000 for 8 = 89°, 60°, 45° and 
25°. The results are also compared with a classical Meyer-Peter&Miiller model. A more 
detailed comparison of the final bottom can be seen in Figure 4(a) 

Remark that for 8 —> 90° we have k —>• 0. Thus, the model reduces to a classical Meyer- 
Peter&Miiller model for values of 8 near 90° as it can be observed in Figure 3(c). For 
smaller values of 8, the effects of the gradient d x (rh\ + rj) play a relevent role in the stress 
tensor. This is shown, for instance, in Figures 4(b) and 4(c), where the value r — r e // 
displayed. The steep profile of d x {rh\ + rj) makes that this term plays an esentiall role in 
the evolution of the test. The smaller the value 8, the bigger the influence of the gradient 
on this test. This is specially true for small values of t. This results in a smoothing of the 
dune profile, and in particular in the advancing front of the dune. 

3.2 Test 2 

The purpose of this test is to study the influence on the definition of the effective shear 
stress in terms of V x ri instead of V^r/ii +rj)/(l — r ). Following Remark 2.1, both terms 
can be related by 


V x i) = -^-Vxirhx + rj) - + h x + h 2 ). (47) 

This means that both definitions coincide only for constant free surface. We recall that 
the model proposed in [26] defines the effective shear in terms of \7 x r]. This model was 
also studied in [40]. 

We set as initial condition 
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t = 0.000 









— Free Surface 




—— Bottom 
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0 2 4 6 8 10 


(a) Test 1: Initial condition 



(b) Test 1: Surface at time t = 2000 



(c) Test 1: Bottom at time t = 2000 

Figure 3: Test 1: Initial condition and evolution for different values of 5 
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t=2000.000 



(a) Test 1: Bottom comparison (Zoom) 

t = 5.000 



(b) Test 1: r — r e ff at time t = 5 (Zoom) 



3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 


(c) Test 1: r — r e ff at time t = 100 (Zoom) 

Figure 4: Test 1: Influence of 5 
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MO, a;) 


(48) 


0.1 + 0.1 ^1 + cos 2^- 4 )) ’ ^ x e t 0 ' 2, 0-6], 

0.1, otherwise. 

MO, a;) + h 2 (0,x) — 1, gi(0, ®) = 1.4, (49) 

The initial condition is shown in Figure 5. 

We set the boundary condition q(t, 0) = 1.4 and we have used the same parameters 
and computational domain described int Test 1. 6 is fixed to 45°. 


t = 0.000 



Figure 5: Test 2: Initial condition 

Figures 6(a) and 6(b) show the difference observed on the surface and bottom when 
the different terms d x (rh\ + rj) and d x rj are used in the definition of r e ff. Due to the shape 
of the surface and bottom, we get that d x rj > + rj) and d x rj < + rj) on 

the upstream and downstream part of the dune respectively which is shown in 6(c). This 
makes that the influence of the gravitational effects in r e // are stronger when d x rj is used 
and results as a more difussed shape observed in Figure 6(b). 

3.3 Test 3 

In this section we present a comparison with experimental data obtained by the Hydraulic 
Laboratory of Escuela Superior de Ingenieros de Caminos, Canales y Puertos (A Coruna 
University) over a channel of 15 m long and 0.5 m width (for more details see [49] and 
[6]). We compare the experimental data with the numerical simulation corresponding to 
the modified Meyer-Peter&Miiler’s model, defined by (42). 

The experimental test was developed by introducing a sand layer in the central part 
of laboratory channel, and inducing hydrodynamical conditions to erode the sand layer. 
The channel has slope of 0.052 %. Sand layer was situated in interval [4.5 m, 9 m], with 
a thickness of 4.5 cm; being media diameter of the grain equals to 1 mm. As boundary 
conditions, an incoming discharge equal to 0.0285 m 2 /s upstream is imposed. The water 
thickness is 0.129 m downstream. 

We have the experimental measurements of the sediment profile at several points for 
t — 10, t = 40 and t — 120 minutes. Then, for the numerical simulation we have considered 
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t = 200.000 



(a) Test 2: Surface at time t = 200 


t = 200.000 



(b) Test 2: Bottom at time t = 200 

t = 200.000 



(c) Test 2: d x (rh\ + r /) and d x r] at time t, = 200 

Figure 6: Test 2: Comparisons at time t = 200 
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as initial condition for the sediment surface a profile produced by an interpolation of the 
data at t = 10 minutes (see Figure 7). The initial condition for the free surface and the 
discharge can be precomputed by considering this profile of the sediment layer and previous 
boundary conditions. 



Figure 7: Test 3: Sediment surface. Initial condition 

For numerical simulation we have meshed the domain with 250 nodes. The CFL is set 
to 0.9. Sediment porosity is set to 0.4. Friction between fluid and bed is modeled using 
a Manning’s law with coefficient equal to 0.0125 over the fixed bed and 0.0196 over the 
sediment layer. 


t= 40 m. t= 120 m. 



Figure 8: Test 3: Sediment surface evolution at t = 40 min. and t — 120 min. 

Comparisons with the experimental data for t — 40 and t — 120 minutes are presented 
in Figure 8. Two numerical simulations has been done, corresponding to the Coulomb 
friction angle 5 — 10° and <5 = 22°. We can observe that in both cases the averaged 
position of the sediment surface is well reproduced by the numerical simulations. We can 
remark that at t = 120 minutes the numerical simulation reproduces correctly the position 
of the discontinuity in the profile of the sediment bed, located at x ~ 10.7. Nevertheless, 
the model does not produce the small sediment bed at the right of the shock. Although, we 
have tested a large range of models to define the solid transport discharge and any of them 
reproduce this advance of the sediment bed at the right of the shock. Which is probably 
a purely tridimensional effect. 
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4 Conclusions 


The SVE model has been deduced in this work by an asymptotic analysis of the Navier- 
Stokes equations. We show that depending on the considered friction law at the interface 
between the fluid and the sediment we can obtain different definitions of solid transport 
discharge. Results have been presented for the case of linear and quadratic friction laws. 
The stop criteria, which is usually represented as a positive part in terms of a critical 
Shields stress, is deduced by including a Coulomb friction law in the internal interface of 
the sediment layer. 

The proposed models present also some advantages with respect to classical ones: they 
have an associated dissipative energy, they are deduced for arbitrarily sediment sloping 
beds and they depend on the thickness of the moving sediment layer, what implies sediment 
mass conservation in general situations. 

Moreover, from the deduction of the model we obtain a modification on the definition 
of the effective shear stress. We obtain that it must be defined in terms of the gradient 
of total pressure, V x (r/ii + b + ho)- By while, classical models consider the gradient of 
sediment bed, V x (b + h 2 ). Nevertheless, by considering this classical definition we cannot 
obtain a model with a dissipative energy. 

Finally, the deduction of the model and the proof of the existing energy dissipation 
equation provides important information on classical models. It allows us to understand 
why classical SVE models do not have an associated energy. Moreover, in order to have an 
associated dissipative energy we may we add an extra term - a term that is deduced in the 
modelling process - in the definition of the shear stress which appears in the momentum 
equation of the Saint-Venant model. This allows also us to deduce that any classical SVE 
exner model, defined by a solid transport discharge which can be written under the general 
form (2), verifies a dissipative energy equation up to second order. 

In the numerical test we presented the evolution of a dune for several repose angle. 
We can see as the morphodynamical form of the dune is well reproduced. The influence 
of the definition of the effective shear stress as we propose in this paper in comparison 
with the classical one is compared in the second test. We can observe that a difference 
appears between both definitions just in the area over the dune, where the free surface is 
not constant. In the last numerical test we have compared with experimental data for the 
generalization of the Meyer-PeterV-Muller model. We obtain that the averaged position 
of the dune is well captured in comparison with experimental data. And the position of 
the shock is well approximated. 
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Appendix 

A Derivation of the sediment transport models 

In this section we develop the derivation of the sediment transport models from an asymp¬ 
totic analysis of the Navier-Stokes equations. These processes are different for each layer. 
In [24] a similar development is performed coupling also shallow water and Reynolds equa¬ 
tions in a multiscale framework but for a transport of pollutant problem. Thus, we consider 
this work as a reference in the development, which in turn is based on the original works 
[30] and [44], 

The derivation is divided in several steps: first, we write the initial system in non- 
dimensional form; next we make an hydrostatic approximation, we assume a suitable 
asymptotic regime and we average the equations out. 

A.l Dimensionless equations 

To begin with, we write the equations and boundary conditions under dimensionless form. 
We set the dimensionless variables, where we must take into account the different nature 
of the fluids in two layers, so we make it separately. Indeed the main property that we 
want to point out is the different order of the velocities, as we have discussed in Section 
2.1. Since we study a coupled system we relate the two characteristic velocities by the 
aspect ratio to indicate that the sediment layer is slower than the fluid layer. 

We denote by H and L the characteristic height and length respectively. To impose the 
shallow flow condition, we assume that the aspect ratio between the characteristic height 
and length is small, as usual we denote it by e — The characteristic velocities are U 
for the layer 1 and U 2 for the sediment layer, consequently, the characteristic times are 
respectively T — jj and T 2 = ^ for each layer. This hypothesis also affects the definitions 
of the Froude and Reynolds numbers. For the sake of clarity we indicate separately these 
variables. 

We consider the “star” notation for the dimensionless variables. 

General dimensionless variables: 

x = Lx* z = Hz* hi = Hh* b = Hb* 

fric: f — p 2 U 2 inCf* fric = p\U 2 iric* H 
Nondimensionalization for layer 1: 

u\ — Uu\ iv \ = eUw* t — ^-t\ pi = piU 2 p\ 

U L U 

Re i =- Fri = —== Hi = piui 

vi 
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Nondimensionalization for layer 2: 


TT * TJ * , L .* P 2 V 2 U 2 * 

u 2 = U 2 U 2 w 2 = £U 2 W * 2 t = —t 2 p 2 = - —p 2 

U 2 £n 

u 2 l u 2 

Re 2 =- Fr 2 — r - n , p 2 = P2V2P2 T m — £ U 2 

122 VgH 

Since we study the coupled system we must take into account that layer 2 is slower than 
layer 1. This related aspect lead us to search for a relationship between the characteristic 
velocities of the two layers holding this property. Thus, we assume that 

U 2 = e 2 U, 


so consequently, T 2 = Rqj — and we can write the dimensionless variables for the 
sediment layer in terms of e as: 


U2—eUu 2 W2=£ A Uw 2 t — 


2 UL „ 2 U 

Re 2 = £ - Fr 2 — £ 2 —== 

v 2 \J gH 


£ 2 U 2 


P 2 = 


£p2»2U 


H 


-P 2 


p2 = P2^2P*2 Tm = £UT t 


We also define the ratio of densities, 


r = — < 1. 
P 2 


Thus, the equations and the boundary conditions written in dimensionless form read as 
follows (we omit the “star” to simplify the notation): 


d tl ui + div x (ui ® ui) + d z {uiwx) - —— div x (-D x Ui) 

d tl Wi + UiV x Wi + W\d z W\ - — 1 — div x (V x uq) 

iXC\ 

1 1 N O 1 1 02 , 1 O 11 

-d 2 (div x ui) - 2—~^—o z wi + — d z pi = - 


£ 2 Re 1 


£ 2 Re 1 


£ 2 Fr 1 


2 > 


(50) 


^ div x tii + d z w\ — 0 . 


£ 2 Re 2 (d t2 u 2 + div x (u 2 <g> u 2 ) + d z (u 2 w 2 )) 

= — V x p 2 + d z (n 2 d z u 2 ) + 2e 2 div x ( / u 2 D x u 2 ) + £ 2 V x (p 2 d z w 2 y, 


< 


£ i Re 2 {d t2 W 2 + u 2 V x w 2 + w 2 d z w 2 ) 

= ~d z p 2 + £ 2 {2d z (/a 2 d z W2) + d z (p 2 di v x u 2 ) + e 2 di v x (/r 2 V x w 2 )) 


2 Re 2 

6 T>f ; 


(51) 


[ div x ?i 2 + d z w 2 — 0. 
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• Conditions at the free surface: 


d tl ri s + u x ■ V x r/ S = wp, 

\--^—D x ui +p x ) d x rj s + -J—V X W! + —^-^—d z u\ = 0; 

y rC6\ J riCi 6 riC\ 

1 2 

— (e 2 V x Wi + d z ui) V x i) s + --—d z wi -pi = 0 . 

riCi rC6\ 

• Conditions at the interface: 

d tl p + u x ■ V x p = wp 
dt 2 V + U 2 ■ V x p = w 2 ; 


(52) 

(53) 

(54) 

(55) 

(56) 


re 2 |V? 7 | 2 ( 2 -^ r L> x 'Ui - p^ + r—(-2(£ 2 V x u;i + d z u x )V x p + 2d z w x ) - rp x 

= --^—U & \Vx'n\ 2 {2.p 2 D x u 2 - P 2 ) - 2e^ p 2 {e 2 V x w 2 + d z u 2 )V x r 7 + 2 e A p 2 d z w 2 - e 2 p 2 
1162 ' ^ 

(57) 

—( -2V x r)(D x ui - d z wi) + (V^i + 4^«i)(l - £ 2 |V x ?y| 2 ) ) 

^ e i \ £ J (58) 

= -frici/l + e 2 |V x ??| 2 ; 

£ 

£V 2 2— ( -2V x r){D x u 2 - d z w 2 ) + (V x u; 2 + \d z u 2 )(l - £ 2 |V x r/| 2 ) ) 

i?e 2 V £ / (59) 

= r- fric-v/l + £ 2 |V x ?d 2 ; 

£ 

• Conditions at the internal sediment interface: 

d t2 h f + u 2 V x r) f - w 2 = — T m ; (60) 


£ 4 M 2 - 2 - ( -2V x rj(D x u 2 - d z w 2 ) + CV x w 2 + \d z u 2 )( 1 - £ 2 |V X ??| 2 ) ) 

Re * V £ 2 / (61) 

= -fric/\/l + £ 2 | V x r/| 2 ; 

£ 

• Condition at the bottom: 


u 2 V x b + w 2 = 0. 


(62) 



A.2 Layer fij: shallow water. 

In this section we obtain the mass and momentum approximated equations for the fluid 
layer. 


Hydrostatic approximation 

We assume e to be small and keep the terms in the system of order zero and one. Then we 
integrate this system in [rj,r] s \. First, we integrate the equation of the horizontal velocity 
and use conditions (52), (53), (54), (55) and (58) to obtain: 


d, 


'ti 


U\dz ) + diva. 


rvs \ 2 

(u\ ® u\)dz ) — ——diva, 

n ) rtei 


D x u\dz ) + V 


Pidz 


2 1 

+ Pi|*=,V^ + -fric = 0{e 2 ). 

(63) 

Now, to get pi we integrate the vertical velocity equation from z to r/ s and use the divergence 
free condition: 

Pi(z) = Pi\z=r, s - - diVa-Mii^J - ^ (z - t} s ) + 0(s 2 ). (64) 


Asymptotic analysis 

To develop the asymptotic analysis, we assume the following hypotheses on the data: 

A.**, fric — efrico, 

Rc\ 

and we consider the development of the unknowns in terms of e: 

hi = h\ + eh\ + 0(s 2 ), u x — u® + eu\ + 0(s 2 ), pi — p\ + ep\ + 0(e 2 ). 

We also introduce for simplicity the notation rf and rf s to write the main order components 
of the interface and the free surface respectively. 

If we just consider the terms of the principal order (e°), we obtain from first equation in 
(50), (58) and (53) that: 

d 2 u\ =0(e); 
d z u ljz=r/ = 0(e); 
d z u ljz=rts =0(s), 

from where we deduce that u\ does not depend on z at first order, so: 

Ui(t, x, z) = Ui(t, x ). 

This is the usual “motion by slices” property of the shallow flows. Under this hypothesis, 
we shall rewrite the expressions above to obtain the final equation for layer 1 at first order. 
First, we write the mass conservation equation from (50) as: 

d t hi + di v x (h\u\) = 0. (65) 


31 



(66) 


From (64) and taking into account the free surface condition (54), the pressure reads, 
Pi(z) = ^ ( z - Vs) - 2ei'oidhvu l j ) + 0(s 2 ). 

Using these values in equation (63), we obtain: 

d tl (Ji\u\) + div x (hi(ui ® u?)) + \^^x{h\) 2 + + fric 0 = O(e), (67) 

where the friction term frico will be specified later. 


A.3 Layer Q?: Reynolds. 

We derive in this section the first order approximation of the evolution equation for the 
layer 2. Due to the more complex dimensionless considered for this layer, it is necessary to 
keep in the derivation process the terms of order zero and one (the unknowns are denoted 
with tilde). Later we will only consider the terms of principal order. 


We write the momentum equations in (51) up to second order : 

-d z (n 2 d z u 2 ) + V x p 2 = 0(e 2 )-, (68) 

Rp 

d zP2 +£ 2 ^l = 0(E 2 ). (69) 

The mass equation for the sediment layer comes from the integration of the incompre¬ 
ssibility equation. Since we decomposed this layer into the static and the moving sediment 
layers, we make this integration in each of them. Thus, for the static layer we have that 
u 2 (z) — 0 for b < z < rjf, and we use conditions (60) and (62) to get 

dt 2 h f — —T m (70) 

Using conditions (56) and (60) we write the mass conservation for the moving layer: 

u 2 dzj = T m . (71) 

The expression for the velocity u 2 will be obtained from the equation (68), but first we 
need to know the value of the pressure p 2 that will be given by equation (69). 


d t2 h m + div. 



Remark A.l Due to the incompressibility equation, we only need to find an equation for 
u 2 . In fact from the mass conservation equation, we have 

rv 

w 2 (z) = w 2 \ z=v - / div x u 2 dz, 

J Z 

where w 2 \ z=11 is given by the kinematic condition at the interface (56). □ 
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Asymptotic analysis 


We consider the following asymptotic regime in order to keep the gravitational and inter¬ 
granular friction effects: 

u 2 — fric/ = efric/ 0 ; 

Thanks to the definition of the dimensionless variables for the layer 2, we have Re 2 = g2 ^ L , 
so for simplicity we will use the notation 

e 3 ZA> 

Re 2 = —, where v 02 = —= 0(1). 

//02 et 

Notice that from the definitions of Re 2 and Fr 2 we have e 2 || f = so for simplicity we 
introduce 

= R&i _ 1 

0 " Frl u 02 Fr'l ' 

As we did for the layer 1 , we develop each unknown in terms of e and introduce the 
notation: 

hm = h° m + eh} m , u 2 = u° 2 + eu\ , p 2 = p° 2 + ep\. 

Note that the integration at this level involves only the moving layer h m . because the 
velocity is assumed to be zero for b < z < rjf. 

So to obtain the pressure we integrate the equation (69) for z G [rjf, rj]\ 

M z ) = P2\z= v - £ Po( z ~ V°)- ( 73 ) 



We use the interface condition (57) to get the value of the pressure at the interface up to 
second order 

P2\ z=v = £—{pi\ z =r, + 2ezyndiv x 'u < j ) ) + 0{e 2 ). (74) 

^02 

Thanks to (66) we get: 

p 2 (z) = £—^h° 1 -£l3 0 (z-ri 0 ). (75) 

* 4)2 f r l 

Thus 

V x p 2 = £— ~^V x h\ + £/3 0 V x ri 0 , (76) 

* 4)2 Pr 1 

that does not depend on z. Taking into account the definition of /3 0 , we can write this 
equation as follows: 

V x p 2 = £—-^{rV x h° 1 + W x r)°). (77) 

* 4)2 Frf 


Next in order to get an expression for the velocity u 2 , we integrate (68) from rjf to 2 to 
find 


p 2 d z u 2 = (p 2 d z u 2 )\ z=Vf + (z- rjf) V x p 2 . 
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We get the value at the internal sediment interface from the friction condition (61) at 
second order: 

vo 2 (H 2 d z u 2 )\ z=Vf = efric /o + 0(e 2 ). (78) 

Then we obtain 

d z u 2 = -— fric / 0 + V xP 2 ~—— • (79) 

h2 ^02 H2 

We integrate again but now from z to r/ to get u 2 . We must take into account that /j 2 is 
not constant, so 


£ Z* 2 1 /*^ £ — Tj f 

U 2 (z ) = u~ 2 | 2=1? H-fric/ 0 / — dz + V a # 2 / -- -dz. (80) 

^02 J v M2 J v M2 

Since we are interested in the principal order approximation, we neglect the two last terms 
—remember that V x p 2 ~ 0(e) from equation (77) — . So we approximate 

u 2 (z) = u~ 2 | z=T) + 0(e). 

If we use the boundary condition at the interface (59), 

(/r 2 5 2 u 2 )p=^ = e—fric 0 + 0(e 2 ) (81) 

^02 

and the previous expression for n 2 d z u 2 evaluated in z — rj we get a relation between the 
friction forces: 

r fric 0 = fric /o + -^(rV^h? + V x r] 0 ) (82) 

where we have used (77). 

At this moment, in order to find an expression for u- 2 \ z=v , and then for u 2 , we must explicit 
the friction terms. 

For the friction at the level z — r]f, we consider a Coulomb friction law, 


fric/ = — (sgn(u 2 ) tan J((<Ti - a 2 )N f ) ■ Nf) (83) 

with Nf the unitary normal vector to the interface z — r]f and S the intergranular Coulomb 
friction angle. To be consistent with the development done before, the asymptotic assump¬ 
tion must be 

tan 5 — etan<5 0 - 

Developing this expression and using ( 66 ) we have: 


r — 1 

fric / 0 = — _ ^Y 2 _ sgn(u 2 ) tan 


(84) 


We consider two possible friction laws for the level z — rj and we derive the corresponding 
models. 
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• Linear friction law at z — r] 

We use a generalized law based in the work [42] (see Remark A.2), that reads 

fric = C(ui - u 2 )\ z=r , (85) 

Remark A.2 The friction law introduced in (85) is based on the work performed in [42]. 
If we take into account the asymptotic regime considered for the viscosities, the coefficient 
C can be written (up to second order) as: 

C — kh m . 


We must also take care of the adimensionalization for this friction term. Thus we assume 
the following dimension and asymptotic to the coefficient C: 

C = UC *; C* = eC°. 

Then we have 

fric 0 = C° (u° - s 2 u 2 \ z=v ) ( 86 ) 

From the friction at z — p ( 86 ) and using (82) we get the value of u 2 \ z=v 

_ 1 o 1 f . 

u 2\z=ri £ 2 e 2(J0^ 11C ° 


f2 U °1 - ( fric /o + + VxV 0 )) 


Finally, an expression for u 2 is deduced: 


u 2 (z) = ^ u° - r -J-qhm ((1 - r)sgn(u 2 ) tan 5 0 + (rV x h\ + V x p°) ) (87) 

Now we use the relation (71) to write an equation for h rn up to first order: 
dt 2 h Q m + div x (^hP m u\- ^ c0 ((! ~ r) s g n M tan £ 0 + (rV x h\ + V x r?°)) ^ =T m . 


dt 2 h° m + div x {^—h 0 ^- r £2C o ~02~ (( x ~ r)sgn(M 2 ) tan£ 0 + ( rV x h\ + V x rj°)j j = T m . 

( 88 ) 

Finally we use (82) to complete the equation for the momentum of the layer 1, (67) 
that reads: 

d tl {h\u\) + di w x (hl(u\ ® u\)) + \^V x {hf) 2 + -^hy xV ° 

1 1 _ f 1 Fri v (89) 

+- jrqh m [( 1 - r)sgn(u 2 )tan5 0 + (rV x h? + V x r]°)J = 0(e). 

The final model is given by equations (65), (70), ( 88 ) and (89). The dimension form of the 
system is given by (17)-(19). 


35 



• Quadratic friction law at z = rj 

We consider now a quadratic law as follows 

fric = CiKui - u 2 )\ z=v \{u 1 - u 2 )\ z =r) (90) 

with C\ an adimensional constant that we suppose of order e, that is, Ci — eC®. Then, 

fric 0 = Ci |u? - e 2 u 2 \ z=ri \ (u° - £ 2 u 2 \ z=ri ). 

We can find a solution for this equation, given by: 

1 


til £ ^ 2 \ z = t ) 




|fric 0 | 1/2 sgn(fric 0 ) 


So we write the value of the velocity u 2 at the interface as follows: 


— Jl o _ 

u 2 \ z=v — 2 u i 


1 


| fric 0 1 1/2 sgn (fric 0 ) 


Thus, (87) turns into 


h 


1/2 


U2 ' [Z £ 2 Ul £ 2 


\V\ 1/2 sgn(V) 


with 


V = (v x (rhi + h 2 + b) + (1 - r)sgn(u 2 ) tan <5^. 


(91) 


Since the relation (82) is valid for any friction laws, the completion of the equation for the 
layer 1 does not change and is given by (89). The dimension form for the complete system 
is given by (17)-(20). 


B Energy balance 

B.l Proof of Theorem 2.1 

In this appendix we prove the result of the energy balance associated to the proposed 
system that is held in Theorem 2.1. Concretely, we prove that the deduced model verifies 
an exact entropy dissipation energy. 

For convenience, we remind here the proposed model given by equation (17) 

( d t hi + div^i = 0, 


d t q\ + div x (hi(ui <8 > uQ) + ^gV x h 2 + ghiV x (b + h 2 ) + = 0, 

< 

d t h 2 + div* ( h m v b yj{l/r - 1 )gd s '\ = 0, 


[ d t h f = —T m . 
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with 


V = (y x {rh\ + h 2 + b) + (1 - r)sgn(« 2 ) tan^ 


and 


for the model ( LF ), 


Vb — v[ LF1 — . u\ — — — V. 

b ^/(l/r-l)gd s l~r 


for the model ( QF ), 


( qf) 

Vb = vf = 


v / ( 1 / r -1 )gd. 

In the following lines, we give the proof of the Theorem 2.1. 


Ul _ \V\ 1/2 sgn(V). 


First we multiply the momentum equation for the layer 1 by ru\ and we use the mass 
equation to obtain: 

^ + rgh\UiS7 x{b + h 2 ) + ghmU 1 'P j = 0 

(a) (b) 


-d t {ghl + hx\ui\ 2 ) + rdiv„ 


h\U\ 


M 3 


+ ghi 


Now we multiply the equation for the layer 2 by V — g{{rhi+h 2 +b)+x (1—r)sgn(-u 2 ) tan 5). 
Note that V X V — gV. 

• For the linear friction law case, 


7 N9f/i2 + Vdw x (h m ui)' - Vdiv x 

(c) (d) "- 


Vi 1 / 7 ' ~ 

' N/ ' 

(e L ) 



= 0 ; 


• for the quadratic friction law case 

£^^ + £div^(/vpii2-Pdiv 2: | h m \j Y^\'P\ 1/2 sgn('P) - l)gd s J = 0. 

(C) ' (d) ' -2---^- L 

(e«) 


We use the definition of V to decompose the term (c) in two parts, (c) = (c)i + (c) 2 : 

(c)i = g(rhi + h 2 + b)d t h 2 

(c) 2 = gx (1 - r)sgn(u 2 ) tan 6d t h 2 — d t (gh 2 x( 1 - r)sgn(u 2 ) tan 5) 
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Now, after some simple calculations, the terms (a) and (c)i gives: 

(a) + (c)i = \gd t ({b + h 2 ) 2 ) + dvv x {rghi(b + h 2 )u 1 ) + rgd t (hi(b + h 2 )) 

Finally 

(6) + (d) = div x (h m u 1 V) 

and the last terms read 

(■ e L ) = —div x (hmy^PV \/(l jr - l)gd^j + ghmj^ \/( 1 / r “ \P \ 2 

or 

(e Q ) = —div x Y~^l^| 1/2s g n (^) \/(!/r- l)s<j Vl 1 / 7, _ l )9 d s |T| 3/2 

So, finally we obtain (43) and (44), where the right hand side in both cases are non-positive. 


B.2 Proof of Theorem 2.2 


In this section we analyze the energy associated to a general Saint-Venant-Exner model 
for a modified friction in the momentum equation of the Saint-Venant system: 


r/pi = gipih^uiluil + -£ m Vgy{hi)K. 

This model is given in equations (l)-(3) that we remind next: 
f d t hi + div x Qi = 0, 

dm + div x + \g h i^ + ghidiv x (h 2 + b) + r/pi = 0, 

( d t h 2 + di v x q b = 0, 


where 


Qb 


Q 1 - <p 

Q = d B \Jg(\/r - 1 )d s , 6 = 


sgn(r) ki 6 m 1 (9 - k 2 0 c )™> (Vd - hy/F c )™ 3 , 

\A d l 


and t — pigip(hi)ui\ui\. 


g{p2 ~ Pi) d l 

To prove Theorem 2.2 we follow the same development above. First, we write the discharge 
qb in the following way: 

^ = CmSgn(r) (Ve - \J~0 C + sgn(r)c/V x (rhi + h 2 + b)j , 
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where 


S m 


1 - ip 


h e m1 (o - k 2 e c )™ 2 (Vo - hVV c )T 


(VO - W + sgn(r )gS7 x (rh\ + h 2 + b )) 


Now, from the definition of 0 and t we write: 


e = so ^ = 

(P2 - Pl)d s ’ \/(l/ r “ 1K> 

with 

TZ = sgn(r) - gV x (rh i + h 2 + b). 

So, finally 


(ft = f m ?ftv / 5 , ( 1 / r - l)d s with v b 


V9<P{h i ) 

V 7 5 , ( 1 / r - IK 


«i — TZ. 


And we can write the equation for h 2 as follows: 


d t h 2 + div x {imVWWui ~ £m'R'Vg(t/r - l)d a ) = 0. 

Note that we have rewritten the evolution equation for the sediment in a different manner 
but no modification has been introduced in it. As we mentioned above, the only modifi¬ 
cation needed to obtain the dissipative energy balance is taken into account in the friction 
term in the momentum conservation equation of the Saint-Venant system. 

First, we multiply the momentum equation of the first layer by ru\ and we use the mass 
equation to obtain: 

'-dt(ghl + hi\u x \ 2 ) + rdiv x (Vtq + rghiUxV x (b + h 2 ) 

(a) 

+ gip(hi)\ui\ A + ZmV gv(hi)Hui = 0 

"-v-' 

(b) 


Now we multiply the equation for layer 2 by TZ — sgn(r)(g(rh 1 + h 2 + b) + xVV), where 
V X TZ — 1Z. Then, 

TZdth 2 + TZdiv x (%mV9V(hi)u^j - 1Zdiv x {^ m TZ y/(l/r - l)gd.^j = 0. 


In the same manner than before, we obtain that 

(a) + (c) = ^gd t ((b + h 2 ) 2 ) + d\\ x (rghi(b + h 2 )u x ) + rgd t (hi(b + h 2 )) - d t (h 2 sg\\(r)\/0 c ) 
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Finally 

( b ) + (d) = div^^UivV^i)^) 

and the last term reads 

(e) = —div x (imp'll \J (1/r - l)c/d s ) + |^| 2 

Then, (45) is satisfied. 
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